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Boltzmann differential operator 
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direct knockout redistribution, MeV~ l 
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High Energy Transport Code 
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High Charge and Energy Transport Code 
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step size for numerical integration 

j(k) r(k) j(k) 

el ' (/ 

integral operators 


integrals from E; to £/+i of terms on right-hand side of 
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LAHET 

Los Alamos High Energy Transport. Code 
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Monte Carlo A T - Particle Transport Code 
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stopping power jth particle 
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solar particle event 

xj. yj ,x 
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Abstract 


A neutron transport algorithm including both elastic and non elas- 
tic. particle interaction processes for use in space radiation protection 
for arbitrary shield material is developed. The algorithm is based 
upon a multiple energy grouping and analysis of the straight-ahead 
Boltzmann equation by using a mean value theorem for integrals. 
The algorithm is then coupled to the Langley HZETRN code through 
a bidirectional neutron evaporation source term. Evaluation of the 
neutron flue nee generated by the solar particle event of February 23, 
1956 , for an aluminum-water shield-target configuration is then com- 
pared with MCNPX and LAHET Monte Carlo calculations for the 
same shield-target configuration . With the Monte Carlo calculation 
as a benchmark , the algorithm developed in this paper sho wed a great 
improvement in results over the unmodified HZETRN solution. In 
addition r a high-energy bidirectional neutron source based on a for- 
mula by Ranft showed even further improvement of the flue nee re- 
sults over previous results near the front of the water target where 
diffusion out the front surface is important. Effects of improved in- 
teraction cross sections are modest compared with the addition of the 
h i gh-e ne rgy b i direc tional sou rve t e rms. 


Introduction 

This paper presents an improved algorithm for the analysis of the transport of secondary neu- 
trons arising in space radiation protection studies. The design and simulation of the operat ional 
processes in space radiation shielding and protection require highly efficient computational pro- 
cedures to adequately characterize time-dependent environments and time-dependent geometric 
factors and to address shield evaluation issues in a multidisciplinary integrated engineering design 
environment. One example is the recent study of the biological response in exposures to space 
solar particle events (SPEs) in which the changing quality of the radiation fields at specific tissue 
sites are followed over 50 hours of satellite data to evaluate time-dependent factors in biological 
response of the hematopoietic system (ref. 1). Similarly, the study of cellular repair dependent 
effects on the neoplastic cell transformation of a C3H10T^ mouse embryo cell culture popula- 
tion in low Earth orbit, where trapped radiations and galactic cosmic rays vary continuously in 
intensity and spectral content about the orbital path (ref. 2), requires computationally efficient 
codes to simulate time- dependent boundary conditions around the orbital path. But even in 
a steady-state environment which is homogeneous and isotropic, the radiation fields within a 
spacecraft have large spatial gradients and highly anisotropic factors so that the mapping of 
the radiation fields within the astronauts tissues depends on the astronaut timeline of location 
and orientation within the spacecraft interior where large differences in exposure patterns that 
depend on the activity of the astronaut have been found (ref. 3). Obvious cases exist where rapid 
evaluation of exposure fields of specific tissues is required to describe the effects of variations in 
the time-dependent exterior environment or changing geometric arrangement . A recent study 
of the time-dependent response factors for 50 hours of exposure to the SPE of August 4, 1972, 
required 18 CPU hours on a VAX 4000/500 computer by using the nucleon light ion section of 
the deterministic high charge and energy transport code HZETRN (ref. 4). In comparison, it is 
estimated that the related calculation with a standard Monte Carlo code such as HETC (ref. 5) 
or LAHET (ref. 6), which are restricted to only neutrons, protons, pions, and alphas, would have 
required approximately 2 years of computer time on a VAX 4000/500 computer. The spacecraft 
design environment, also requires rapid evaluation of the radiation fields to adequately determine 
effects of multi parameter design changes on system performance (refs. 7 and 8). These effects 



are the driving factors in the development and use of deterministic codes and in particular the 
HZETRN code system which simulates 56 naturally occurring atomic ions and neutrons. 

The basic philosophy for the development of the deterministic HZETRN code began with the 
study by Alsmiller et al. (ref. 5) using an early version of HETC, wherein they demonstrated 
that the straight-ahead approximation for broad beam exposures was adequate for evaluation of 
exposure quantities. Wilson and Khandelwal (ref. 9) examined the effects of beam divergence 
on the estimation of exposure in arbitrary convex geometries and demonstrated that the errors 
in the straightahead approximation are proportional to the square of the ratio of the beam 
divergence (lateral spread) to the radius of curvature of the shield material. This ratio is small 
in typical space applications. From a shielding perspective, the straight-ahead approximation 
overestimates the transmitted flux and the error is found to be small in space radiation exposure 
quantities. Langley Research Center’s first implementation of a numerical procedure was 
performed by Wilson and Lamkin (ref. 10) as a numerical iterative procedure of the charged 
components perturbation series expansion of the Boltzmann transport equation and showed 
good agreement with Monte Carlo calculations for modest penetrations to where neutrons play 
an important role. The neutron component was added by Lamkin (ref. 11) and closed the gap 
between the deterministic code and the Monte Carlo code. The resulting code was fast compared 
with the Monte Carlo codes but still lacked efficiency in generating and operating with large 
data arrays which would be solved in the next generation of codes. 

The transport of high-energy ions is well adapted to the straightahead approximation. In 
fact, the usual assumption that secondary ion fragments are produced with the same velocity as 
the primary' initial ion (ref. 12) is less accurate than the straight-ahead approximation contrary 
to intuition (ref. 13). The Boltzmann transport, equation for the particle fields <pj(x,E) is given 
for the straight-ahead and continuous slowing down approximations as 

37 - ^(£) + ^(£)| 4 >j(*,E) = £ <rjk(E£') h(^E') dE' ( 1 ) 

L J k JE 

where x is the depth of penetration, E is the particle kinetic energy, SA E) is the particle stopping 
power, (Tj(E) is the macroscopic interaction cross section, and <7 j^{E,E > ) is the macroscopic 
cross section for particle j of energy E produced as a result of the interaction w’ith particle k 
of energy E 1 . It has been customary in codes developed at Langley to invert the differential 
operator and implement it exactly as a marching procedure (ref. 14). The remaining issue has 
been to approximate the integral term on the right-hand side of equation (1). The formulation 
of the code to approximate heavy fragments was facilitated by assuming that their fragment 
velocity is identical to that of the primary ion velocity. This assumption is inadequate for the 
description of the coupled nucleonic and light ion components. A computationally compatible 
nucleonic transport procedure was developed by Wilson et al. (ref. 15) and agreed well with 
exposure quantities evaluated by Monte Carlo transport procedures (ref. 16). The transport 
of the nucleonic component was developed by assuming the midpoint energy, within the step 
size h, was the appropriate energy to evaluate the integral term. Thus, the residual range of 
the proton will reduce by hf 2 before the interaction, and the secondary proton residual range 
will reduce by h/2 before arriving at the next marching step. Neutrons show no loss in residual 
range as their stopping power is zero. This choice was shown to minimize the second-order 
corrections to the marching procedure (ref. 17). Although reasonable agreement on exposure 
quantities from Monte Carlo calculations was obtained, the resultant neutron flux at the lowest 
energies was substantially below the Monte Carlo result in the range of 0.01 to several MeV 
and required improvement (ref. 18). Analysis revealed that, the problem was in the rescattering 
terms in which the number of elastic scattered neutrons was underestimated numerically, which 
must be addressed as suggested by Shinn et al. (ref. 18). 


2 



The issue of evaluation of the integral term of the Boltzmann equation for the elastic 
scattering was developed in a prior report (ref. 19). In reference 19, a multiple energy 
group (multigroup) method based upon a mean value approximation to the integral terms for 
transporting evaporation neutrons showed vast improvement in the low-energy neutron spect ra. 
Now that the elastic scattering events are adequately represented, consider the addition of 
improved estimates of nonelastic processes on the neutron transport solution. In the present 
paper, nonelastic processes are added to the algorithm developed by Heinbockel, Clowdsley, and 
Wilson (ref. 19). The code is then modified to account for the high-energy neutron production 
at backward angles. Improved neutron interaction cross sections with less dramatic changes on 
the neutron spectrum are introduced. 

Formulation of Transport Equations 


Define the linear differential operator as 

B\6) = 


h-ksm+'t* 


dx dE 
d<f>{x.E) d 


<p(x,E : 


- dx dE l S A E ) + a A E ) *(*.£) 

and consider the following one-dimensional Boltzmann equation from reference 20 


B[<t>j] = Y, r Vjk&tf) dE’ 


( 2 ) 


(3) 


where <f>j is the differentia] flux spectrum for the type j particles, Sj(E) is the stopping power 
of the type j particles, and crj(E) is the total macroscopic cross section. The term 
a macroscopic differentia] energy cross section for redistribution of particle type and energy, is 
written as 

a jk (E,E') = Pd °A E ') fjk,iAE,E') (4) 

li 

where / ^(E <E l ) is the spectral redistribution, <7$ is a microscopic cross section, and pj is the 
number density of [3 type atoms per unit mass of material. The spectral terms are expressed as 


- fll.n + fjk,0 + 


ijk,d - J jk.d ■ 


(5) 


where /jj. j represents the elastic redistribution in energy, ^ represents evaporation terms, 

and fjfc j represents direct knockout terms. The elastic term is generally limited to a small 
energy range near that of the primary particle. The evaporation process dominates over the 
low-energy range (E < 25 McV), and the direct cascading effect dominates over the high-energy 
range (E > 25 MeV), as illustrated in figure 1 with data from reference 20. 

Equation (3) is then written for j = n as 

Bfc.1 = E T E n -v<£'> (/'.La + r„ kJ + rtj) m*j?) («) 

k JE d 
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which is expanded to the form 


%J = r T.'* 'If* 1 ' Ots + Inn:! + HU) M*.2) *2 

+ E /, E Pa °a(£’) (/Xa + + /«m) ^ (7) 

kfr> Jt ji 

Define the integral operators as 

/ oc> 

E pa °a(E') f«l 3 <f>(x,E r ) dE ' (8) 

J 

4 *V] = f Y.Pa <ra(E') fukjj <t>(x,E') dE' ( 9 ) 

J E 0 

I { d k) y>] = r E Pa <ra(E') fit .J 4>(x,E’) dE' ( 10) 

j e :j 

where k = n denotes coupling to neutron collisions, k — p denotes the neutron source from 
proton collisions, and similarly for other ions. When considering only neutrons and protons, 
equation (7) can be written in the linear operator form as 

B[<pn] = il | * [<pn] + l{ \<pn] + Ij ^ [<j>n] + 4j ^ [<j>p] + 1^ \<Pp] + I^J \<f>p] (11) 

Note that does not contribute to the neutron field because protons cannot produce 

neutrons through elastic scattering and therefore equation ( 11), with</>„ replaced by <f>, is written 
as 

BWi — ^e) ^ + Id V] + + l\( \<t>p] (12) 

In reference 19, we assumed the evaporation source to be isotropic and evaluated the t ransport in 
forward and backward directions by using the straight-ahead approximation for elastic scattering 
and found improved agreement with Monte Carlo calculations. The first step in this study Is to 
add effects of nonelastic events into the transport process. 

Assume a solution to equation (12) of the form <j> = <j> t + <^ ( y, where <j> t is the solution for 
evaporation sources and contributes over the low-energy range and (pj is the solution for the 
direct knockout sources and contributes mainly over the high-energy range as suggested by 
figure 1. Substitute this assumed solution into equation (12) and find 

B[^) + B(<t> d ] = /<;%] + /Wferf] + /j B) [*e] + l} n) [0 d ] 

+ I/ \<t>t] + \<j> d ] + p p \<j>p] + 1 [o p ] (13) 

In reference 19, the terms l\ n \<j> t \ and Ij " ^ [<i>, ] were set to zero to consider only elastic scattering. 

This allows estimates of the elastic scattering effects on the transport of evaporation neutrons. 
In contrast, these terms are retained and they demonstrate nonelastic effects on the transport 
of evaporated neutrons. This change also allows flexibility in further improving the HZETR.N 
code as shown later in this report.. As in reference 19, we assume that <j> d was calculated by the 
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HZETRN program because for the direct cascade neutrons the straight-ahead approximation is 
valid. Consequently, <p c { is a solution of the equation 

B[M = 7<rW] + Id%d] + l { j\op] ( 14) 

This assumption simplifies equation (13) to the form 


B[*A = fffa] + ll U) [4>e] + I { d U) m + + li P %] 

The elastic and reaction scattering terms are defined as 

*r,!i = Pd *rp(E') \r nk JE.E') + fitJE.Ef)] 


(15) 


with units of cm^/g-MeV, and for neutrons the stopping power Sj(E) is assumed zero, and 
therefore, equation (15) reduces to the int.egro-differential transport equation with source term 
as 

= V r [<r*ji(E,E!) + <r r j(E,E!)] <p,{x,E?) dEf + g(E,x) (1(5) 
^ J E 


d_ 

dx 


+ a(E) 


Equation (16) represents the steady state evaporation neutron fluence d,(x,E) at depth x and 
energy E . The various terms in equation ( 16) are energy E with units of MeV, depth in medium 
is x with units of g/cm 2 , <f> t \x.E) (in partieles/em 2 -MeV) is the evaporation neutron fluence, 

and g( E, x) — + Ie[<j>p\ (in part icles/g- MeV) Ls a volume source term to he evaluated by 

the HZETRN algorithm. Equation (16) is further reduced by considering the neutron energies 
before and after an elastic collision. The neutron energy E after an elastic collision with a 
nucleus of mass number At.., initially at rest, is from reference 21 


E= E' 


+ 2.4j^cos0 -f 1 
( A T;i + l) 


(17) 


where E f Ls the neutron energy before the collision, Aj> is the atomic weight of t he 3 th type of 
atom being bombarded, and 0 is the angle of scatter. Note that for forward scattering 0=0, 
E — E f , and for backward scattering 0 = tt, E — E r a j. where is the ratio 


<*& = 


At.* + 1 


(18) 


which Ls a constant less than 1. Therefore, change the limits of integration for the elastic 
scattering term in equation (16) to [E,E/& 3 \, which represents the kineticallv allowed energies 
for the scattered neutron to result in an energy E. Equation (16) then is written as 


aT + ^> 


<p.{x,E ) 


Y. [ E/01 <r>AE,E , Un{x.E , )dE' 

3 JE 

+ E r ^rA E ^) 4>A*rf)dE? + g{E<x) 

3 


( 19 ) 
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The quantity <r(E) has units of cm 2 /g and is a macroscopic cross section given by 

a(E) = o$(E) v'^E) (20) 

3 3 

where pj is t he number of atoms per gram, o^(E) is a microscopic elastic cross section in units of 
cirr/atom. and a^( E ) is the corresponding reaction cross section. Other units for equation (16) 
are obtained from the previous units by using the scale factor representing the density of the 
material in units of g/cnA 

Mean Value Theorem 

Throughout the remaining discussions the following mean value theorem is used for integrals. 

Mean Value Theorem: For <ft(x,E) and f(E ) continuous over an interval a < E < b 
such that ( 1 ) <j)(x,E) does not change sign over the interval (a, 6 ), ( 2 ) is integrable 

over the interval (a, b), and (3) f(E) is bounded over the interval (a ., 6 ), then there exists 
at least one point c such that 

f(E) 4>(x,E)dE= /(e) f f> (p(x,E) dE (a<<<b) (21) 

Ja Ja 


In particle transport this mean value approach is not commonly used. In reactor neutron 
calculations, an assumed spectral dependence for <p(x,E) is used to approximate the integral 
over energy groups. The present use of the mean value theorem is free of this assumption; thus, 
more flexibility is allowed in the HZETRN code, and the result is a fast and efficient algorithm 
for low-energy neutron analysis. 


Multigroup Method 

To solve equation (19), partition the energy domain into a set of energies {Eq,E\,. 

Ely £ 7 + 1 , . . Consider first the case where there is only one value of /? which represents neutron 
penetration into a single element material and let <p t be denoted by <f>. Equation (19) is integrated 
from Ej to £/_ j_i with respect to the energy E to obtain 



d<j>{x,E) 

dx 



a{E) <j>(x,E) dE = / 6> / + I ri + £, : 


where 


The quantity 


[ E i + 1 [E/<* 

I s j = / / <r»{E,Ef) Mx,Ef) dE* dE 

JEj JE 

f E i+\ f 30 

/,. i = / c r (E,E ) 4>(x,E') dE' dE 

JEj JE 


r E i+\ 

;i= / g(E,x 

Je, 


) dE 


4> 


f E i + 1 

;(*) = / <ttx,E) dE 
Je, 


' Ei 


( 22 ) 


(23) 

(24) 


(25) 


(26) 
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is associated with the ith energy group (Ej y E ;+\ ), so that ^ — *_£• $j(x) represents an average 

fluence for the ith energy group. Then equation (22) can he written as an ordinary differential 
equation in terms of <J>,(x) as follows. In the first term in equation (22), interchange the order 
of integration and differentiation to obtain 

[ Ei + 1 mx£) dE _ d*£x) (27) 

If. dx dx 


By using the previously stated mean value theorem for integrals, 
can he expressed as 



o' <j)(x,E) dE = W^j(x) 


the second term in equation (22) 


(28) 


where <f = rr[E; + ft(E; + \ — E;)] is a mean value associated with some valueof 6 between 0 and 1. 


For the term l s j in equation (23), interchange the order of integration as illustrated in 
figure 2. The integration of equation (23) depends upon the energy partition selected. For 
example, figure 2(b) illustrates an energy partition where Ej . yi < Ej/ o : for this case we can 
write equation (23) as 


I,.i = dE dE + / / 

J F?=Ej JE=Ej J E'—Ej + \ d E—Ej 

f Ej+ 1 / « fE f+1 

+ / / //, dE dE' 

JE' = EJ< * J E=a E' 


//, dE dE' 


(29) 


where H s — cr s ( E .E 1 ) <s(x. E ' ). Figure 2(c ) depicts the case where Ej+ j — Ej/n exactly for all i. 
In th is special case, equation (23) reduces to 


/•£*+ 1 r£ , r E i+ 1/« r^+i , 

/ ■= / / H s dEdE'+ / H,d.EdE' 

JE l =E t JE=E i JE' = E i+] JE=aE' 


30) 


The selection of an energy partition can lead to two or more distinct energy groups associated 
with each interchange in the order of integration. For example, see figure 3(a). 

The evaluation of equation (24) is somewhat more complicated. As an approximation, we 
assume there Ls an energy E . y such that <f>{x. E) can be taken as zero for all E > E y. In this 
case, equation (24) can be written as 


/•£;+i f E _ — r f^jE 1 f^iE l , 

= L r L r »<■ dE dE + L L L H ‘ iE iE 

J E' = Ej J E= Ej j_ ( + , JEj J Ej 

where H r = a r {E,E') <p{x,E'). For example see figure 3(b). 

Equations (30) and (31) may then be written for the case where E, + \ = E, / o as 

hi = [ E ' <? S (E,E*) dE 4>;(x) dE+ <r,(E,E* +l ) dE * i+l (x) dE 

JE; J<\E* . 


rE j+] rEf n 


and 


f E* -hh f Ej + 1 

/,.;=/ M £.£?)<*£ $*(*)<*£ + zZ / Vr(E,Ej) dE 4>j(x)dE 

J E, j=,-+l ^ 


(31) 


(32) 


(33) 
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where E* = 9\(Ei+\—Ei) and Ei+\ = 62 ( 8 i+ 2 ~ £V+i ) for some 6 \ and 0% such that 0 < 61,62 < 1 
by once again using the previously stated mean value theorem. The special partitioning of the 
energy as illustrated in figure 2(c) enables us to obtain from equation (22) a system of ordinary 
differential equations of the form 



*0 


'a n 

a V2 

«13 • 

• a 1N ' 


*0 


' £0 1 

d 

<*>1 



022 

023 ' 

■ a 2N 


*1 


6 

dx 




- 0 - 

033 ' 

a A \ v . 


$ A _1 _ 

+ 

1 

7 

1 


(34) 


where each equation is associated with an energy group. This is where the term multigroup 
method originates. In equation (34), the coefficient matrix has the elements 


i E ' 

a i,; = I 
Je, 


[<T s (E,Ef) + a r {E,E\)] dE - *(£*) 


r E i + 1 4 r ^>+-2 

«m+ 1 = / *.{E t Et +l )dE+ / <r r {E,E * +] ) dE 

Jo,E *i+\ JE t+] 

f E *+ 1 

«U+j = / <tAE,E* +j ) dE ( j = 2,3,...) 

J E:, ; 


Further assume that for some large value of TV, <F, equals 0 for all i > TV. This assumption 
gives rise to the following system of ordinary differential equat ions: 


dy 

dx 


= Ay + £ 


subject to the initial conditions |7(0) = 0. Here y is the column vector of <!>,' values, 
col($0i $ 1 , . . . , $ v-l)> the matrix A is an TV by N upper triangular matrix, and £ is the column 
vector co1(4q,^j, . . . ,£v-l)- This system can be solved by using back substitution. In a similar 
manner, the integrals in equation (29) and (31) can be evaluated for other kinds of energy par- 
titioning, and a system of equations having the same form of equation (34) obtained. How the 
elements of the matrix A are calculated will depend upon the elastic scattering as determined 
by the type of energy partition. (See, for example, fig. 3(a).) 

For our purposes the system of equations (eq. (34)) is used to discuss some of the problems 
associated with the multigroup method. Of prime concern is how an energy grid is to be 
constructed and how this energy grid controls the size of the matrix in equation (34). Consider 
the construction of the energy partition 


a a'2 qN 

where Eq = 0.1 MeV, for the selected elements of lithium, aluminum, and lead. Table 1 
illustrates integer values of TV necessary to achieve energies greater than 30 MeV. These values 
of TV represent the size of the matrix associated with the number of energy groups. The value 
Eq — 0.1 MeV, in terms of human exposure, represents a lower bound where lower energies are 
not important. The value of 30 MeV represents an upper limit for the evaporation particles and 
could be adjusted for other source terms. 


8 


Tahir 1. Energy Partition Size N 


Element 

o 

N 

0.1/q' v 

Lithium 

0.563 

10 

31.53 

Aluminu m 

0.862 

39 

32.75 

Lead 

0.981 

298 

30.38 


Observe that for energy partitions where Ej+\ < E/jo the values of N are larger, and 
if Ei+ 1 > E}/ a the values of N are smaller. The cases where Ej+\ > Ej/a give rise to 
problems associated with the integration of the elastic scattering terms over the areas and 
,4*2 of figure 2(d) when the order of integration Is interchanged. In this figure, the area A\ is 
associated with the integral defining <!>,, and the area ,4*2 is a remaining area associated with 
an integral which is some fraction of the integral defining $;+i which is outside the range of 
integration; therefore, some approximation must he made to define this fractional part. This type 
of partitioning produces errors, due to any approximations, but it has the advantage of greatly 
reducing the size of the N by N matrix A at the cost of introducing errors into the system of 
equations. A more detailed analysis of the energy partit ion can be found in reference 22. 

The case of neutron penetration into a composite material gives rise to the case where there 
is more than one value of in equation (16). In this special case, equation (23) becomes 


hj = Y,J E f E °*d E ' E ’) n xE '> dE ' dE 

and equation (24) becomes 

j: rx,. <j>{x,E’) dE’ dE 


(35) 


(36) 


In order to avoid the errors introduced when an energy grid is selected such that Ej+i > Ej/a. 
we select a = max(» 1 ,o. 2 ,. . . , Oj) and construct the energy partition where E i+ i = Ej/a so 
that Ej + 1 < E;/a,j for all ft. Obtain a system of differential equations having the same upper 
triangular form but with elastic scattering contributions for off-diagonal elements. Observe that 
for some arbitrary energy grouping we have, for the element hydrogen, a case where the value 
of a.j is zero and Ej/aj is therefore infinite. In this situation, we must integrate over many 
energy groups. In this case, the area of integration is similar to that, shown in figure 3(b). 
For any composite material, depending upon the selected energy partitioning, some type of 
approximations must be made when the order of integration is interchanged in equation (35). 
Also the problem of selecting the mean values associated with each of these integratioas exists 
and now addressed. 


Mean Value Determination 

A realistic test case wassolved analytically and numerically (with and without the multigroup 
approximation) for which the mean values were found empirically for several single element 
materials (ref. 19). The values determined are 

E* = E i + 6 l (E i+l -E i ) 

£ 7+1 = E i+i + e 2 (Ei + 2- Ei+i) 
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where 


and 


where 


( 7! +m n (E -E n ) <$! 
*1 = | Tl + m 12 (£-£ 11 ) -8 X 
{ 73 +m lz {E- £ 22 ) -Si 

(E>E u) 
( E 2 2 < E < En) 
(E < E<2 2 ) 

0-2 = < 

' 72 + oi 2 i(E - E 1 1 ) 
72 4- rn 2 2{E - E\\) 

. 74 + m- 23 (E - £-22 ) 

( E > E\ i ) 
(£22 < E < Eu ) 
(E<E 2 2) 

71 = 0.93 

72 = 0.90 

73 = 0.30 

74 = 0.27 

mu = 0.0030485 
mi 2 = 0.2490258 
mi3 = -0.3937186 
£n = 3.037829 

77J2i = 0.004355 
m 2 2 = 0.249026 
m 2 3 = -0 .255920 
£22 = 0.5079704 


and 5\ is 0.0 for lead, 0.02 for aluminum, and 0.075 for lithium. These values of 6 for the mean 
value theorems were determined by trial and error so that the multigroup curves would have the 
correct shape and agree with the numerical solution. These selections for the mean values are 
not unique. 

Application to Evaporation Source in AI-H2O Shield- Target Configuration 

Apply the previous development to an application of the multigroup method associated 
with an aluminum- water shield-target configuration. In particular, consider the case w'here 
the source term g(E,x), in equation (16), represents evaporation neutrons produced per unit 
mass per MeV and is specified as a numerical array of values corresponding to various shield- 
target thicknesses and energies. The numerical array of values is produced by the radiation code 
HZETRN developed by Wilson et al. (ref. 4). This numerical array of source term values is 
actually given in the form g(Ej ) Xj,yj s .) in units of particles/g-MeV, where yj, represents discrete 
values for various target thicknesses of w r ater in g/cm 2 , xj represents discrete values for various 
shield thicknesses of aluminum, also in units of g/cm 2 , and E; represents discrete energy values 
in units of MeV. These discrete source term values are used in the following w r ay. Consider 
first the solution of equation (16) by the multigroup method for an all-aluminum shield with no 
target material, that is, target thickness = 0. The HZETRN program was run to simulate 
the solar particle event of February 23, 1956, and the source term g(Ei,Xj,y^) associated with 
an aluminum- water shield-target, configuration was generated for these conditions. Using this 
source term, we solved equation (16) by the multigroup method. 

For a single shield material with only one value of /?, equation (16) becomes 

r E/« 

<f>(x,E) = / a g (E,E')<p(x,E’)dE' 

Je 

+ [ <r r (E,E’) <j>(x,E') dE' + g(E,x) (37) 


■Z- + <r(E) 

ox 
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where an integration of equation (37) from FJj to E;+ 1 produces 


f£i+ 1 QA, f£i+ 1 

/ - -dE + <t(E) <j)(x,E) dE 

./£•, J/., 

= / / ax(E,E ) <p(x,E') dE dE 

JEj JE 


+ 


E 

r E i + 1 r 30 . . 

/ / a, •( £■,£/) dtf dE+ / <7(£,*) <*£ 

./£,■ j E: 


( 38 ) 


We define the quantities 


r E >+\ 

4>, - j o(x,E) dE 


r E i + 1 

«■■=/* 9 


(£,*) 


(39) 


and interchange the order of integration of the double integral terms in equation (38). If the 
energy grid is chosen so that Ej + j = Ej /a , a mean value theorem is applied to obtain the result 

(/$ r Ei + ] r E ' i , , 

ZLL + 7r^i= / <r„( £,£■') dE <t>(z,E?) dFf 

dx Je i JE-Ej 

f E i+V f E i + 1 . . 

+ / / cts(E,E ) dE <p{x,E ) dE' 

■I E i+ 1 J E=<y\E' 

+ f ,+l f trAE.Ef) dE<f>(x,Hf) dEf 
JEi Je=E, 

f E J+\ [ E ’+ 1 . ftt /• 

+ >/ / <r r (E,E?) d.E<i>(x,E') dE 1 ’+& 

jtT+l Je j Je > 


(40) 


over the energy group Ej < E' < Ej+\. The first double integral in equation (40) represents 
integration over the lower triangle illustrated in figure 2(c). The second double integral in 
equation (40) represents integration over the upper triangle illustrated in figure ‘2(c). Define 


(41) 


and 



r E ' 

g\(E') = 

I a * 

« 


g-2(E') = 

( EiM , 

- 

' E=a\ E' 


r { (£) = f (t v (E^) dE 
JE=K i 

fEi + 1 , 

r-2i J n(Kn) = J e <rAE,E? m )dE 


( 42 ) 
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then employ another application of a mean value theorem for integrals to write equation (40) in 
the form 


— + cr<*>; = 9\ [ E i + 0\(E,+ \ - Ei)]<$; + g‘ 2 [Ej + \ + 0% (Ej+ 2 - £V+l)]$/+l 

.V-l 

+ n [Ei + 0i(E i+ i - Ej)) *,• + r2ij [ Ej + 0 2 (E j+l - Ej )] ^ (43) 

J=*'+1 


Tliis produces the matrix coefficients associated with the energy group Ej to Ej+ \ so that 


a u = 9i + n ~ ^ 1 

^ m+ 1 3 P2 d - r *2/,/’+l / (j = . . .) (44) 

a M+j ~ r 2/,/+j J 


In this way, the diagonal and off-diagonal elements of the coefficient matrix in equation (37) are 
calculated. 

For a compound target material made up of more than one type of atom, we modify slightly 
the solution technique given in reference 19. For a target material comprised of component 1 
and component 2, there are two values for a. A value a\ is determined for component 1 and a 
value a*2 is determined for component 2 of the compound material. In this case, equation (37) 
takes on the form 


|r + -’<*') 


[E/u\ 

<Kx,E)= j <r S] (E.E 1 ) <j>(x,E? ) dE' 

rE/iX'2 f oo 

+ J a 8 . 2 (E,E r ) <p{x,E’) d.E f + j a ri (E,Ef) dE' 

+ f\r,(E,E') 4>(x,E') dE' + g(E , 

d Fj 


>) 


(45) 


where cr S] and a s ^ are scattering terms and a r] and <J r ^ are reaction terms associated with the 
respective components of the compound material. These terms are calculated in the HZETRN 
code. We consider two cases. Case I requires that the E/a ) line be above the E ja\ line. In 
case ii, a*i equals 0 (the hydrogen case), and the limit of integration for the second integral goes 
to infinity. Each case is considered separately. 

For case I, we assume that a\ > a*j > 0 and select the energy spacing Ej+ \ = E;/a\. We 
then proceed as we did using the single component shield material. Integrate equation (45) 
from Ej to Ej+\ and interchange the order of integration on the double integral terms. Define 

£/ — $(E,x) dE and obtain the equation 


n<P; 

+ = H pi + H 12 + H-2i + Hn + An + A'i2 + A21 -f A22 + & (46) 

where H& \ and Hg 2 represent the elastic scattering caused by collisions with 0 type atoms and 
E f j\ and represent the nonelastic scattering. Note that Hgi and Kg\ are integrals over the 
energy range (Ej, Ej+i) for 0 — 1, 2, and H ^ an( J K $2 represent integrals over higher energies. 
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The integrals H\\, H\ 2, K n, and A'12 are the easiest to evaluate because of the exact spacing 
of the energy partition. These integrals have the forms 


f E i + 1 f e . 

H U = / <rs](E,E ) dE 

JEj J E=Ej 


7 <p{x,E') dE' 


[ E i+'l f E i + 1 . . , 

H\2 = / / a >S] {E,E')dE <p{x.E?)dE? 

■>E i+ x JE=,x 1 E' 


/■£■<+ 1 /■£' , . 

A’ll = / / (Tr,{E,E ) dE o(x,E ) dE 

>E, JE=Ei 

N ~ l r E j+ , ,£, + l 

A- 12= 5 ] / / VrxiE.E') dE <f)(x,E') dE 1 

* ■ 1 1 v f 1 


j=/+l 17 

Here the first subscript represents the material component. A second subscript of 1 represents 
integration over the lower triangle in figure 2(c). A second subscript of 2 represents integration 
over upper triangles, like figure 2(c), or higher rectangles, like figures 3(a) and (b). Defining the 
terms 


*1 (d)( £ ") = 

/ a H AES') dE 

JE=E, 

(tf= 1,2) 

(51) 

h 2( !,(£’') = 

f E i+ 1 . 

/ <t,,i(A\A') dE 

Je=,x I E' 


(52) 

i(,rf )( /^ ) = 

/ <A- E,E') dE 

JE=Ej 

(0 = 1.2) 

(53) 


/ a r j(E<E')dE 

JE-E ; 

(0= 1,2) 

(54) 


and using the mean value theorem for integrals we obtain from equations (47) through (50) 

Hu = /M(i)[£V + 0i(£i+i - 

H\2 — h-2{ 1)[^/+1 + ^2(^+2 ~ £|+l)l^i+l 

A n = k l{l) [E; + B[(E i+ i - Ei)]*; 1 (55) 

.V-l 

/VI 2 = E /-2(1 )[Ej + 0' 2 (E j + ] - Ej)]*j 
j=i + 1 > 

where O\, 0‘2 and define intermediate energy values associated with the mean value theorem. 

The integrals H-2\ and H ? 2 are associated with integration limits {E,E /n 2 ) and energy 
intervals dictated by the selection of a\ for determining the energy spacings. The integral 
H- 2 \ is associated with the triangular area shown in figure 3(a) and takes the form 


H-2i = [ ' +] \ a,.,(E\E')dE p(x,E') dE' 

■IE; Je=E ; ‘ 
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The integral #22 is associated with the remaining shaded area shown in figure 3(a). This 
remaining area is made up of a number of rectangles, trapezoids, and triangles. We approximate 
the integral over each of these rectangles, trapezoids, and triangles as a fraction rjj of the integral 
over the whole rectangle with area A - = (E;_ m — Ei)(Ej+i — Ej ). Integral # 22 , therefore, takes 
the form 

[ e j + 1 r E i + 1 , , , 

H-n = 2_] V / / ^ 2 (£,£') ^ *(*,£) (57) 

j=i+i * /£ j J E > 

where 


.4,/ - 0.5(» 2 £' j >i - E, ){ E J+ \ - Ej/a-j) 

.4 Jj 

_ . 0-5[(Jg i+ i - aiE i+x ) + {E i+ 1 - o 2 ff,)l( ff, + i - Ej) 
Vj - <; 7(77 

0 .5( Ei + i/a-2- Ej )( £V + 1 - a? E , ) 



A IJ 


Defining the term 


/i 3 ( 2 )(£' 


. f E i + 1 

')= ^ 

JE=E; 


( £ >+' < &) 1 
- * 

(i<^ 

(%- < £■>) J 


< ^ < ^+0 
Ej ~ r-1 ^ r~< Ej± \ \ 

< £ ;+i < “srj 


2 {E,Ef)dE 


H' 2 i and H \ 22 can be written as 


(58) 


(59) 


H‘2i — 

N - 1 

7^22= ^3(2)*/ 

J=i+1 


(60) 


Similar to An and A'i2« integrals A 21 and A '22 are given by 


74-21 = 7r 1(2) $, : 

A-l 

/v 22 = X] *2(2)^j 

j=/+l 


(61) 


The coefficients for our system of differential equations (eq. (34)) are then given by 

a i,i = * 1 ( i ) 4- 7ii(2) + * i(i) + 7ri( 2 ) - 7r 1 
°M+1 = *2(1) + Vi+ 1*3(2) + *2(1) + *2(2) 

°M+2 = **'+ 2*3(2) + *2(1) + *2(2) > 

a M+3 = T //+3*3(2) + *2(1) + *2(2) 


( 62 ) 


where evaluation at the appropriate mean energies is implied. 


14 


In case II, the second component is hydrogen which means that r >2 = 0; therefore, one of the 
limits of integration becomes infinite. We once again let oq determine the energy spacing and 
integrate equations (47) through (50) over an energy interval (Ej^Ej- j_i ) which is determined by 
the E f = E/a \ line. Using the definitions given by equations (39) we integrate equation (45) 
over the interval ( E ; , Ei +\ ) and then interchange the order of integration in the resulting double 
integrals to obtain 

™ 

dx 


-d- + — H\ + H>2 + A I + A -I + £>i 


(63) 


/here 


fff = / ?+1 f { (E,E f ) c iE<j>(x.E f ) dE 1 

J Ej JE-Ei 

f E i+‘2 fEi - fl ( , # 

+ / / <r sx {E,E) d,E<p{x,E f ) dE f 

J Ej+ j J E—(.\ \ E f 

H*=[ [ (TkJ E,E') dE <f>(x,E?) dE 1 

“ JEj J Ej 

+ Z] f E+J+] f E + l (T ^( E ' E ''> dE<f>(x,E?) dE’ 

j_^JE j+ j JEj 

' + 1 r & 


(64) 


(65) 


/• £'i+ 1 r c- , 

A'* = J J <r r ,u(EX') dE <p(x,E ) dE’ 


'E f JEj 
JL rE 


+ 


Z] /. ' +;+ ' dE <i>(x,E') dE' 


(66) 


j = 1 ' E ' + ) 


where for all A r * greater than some integer N > 0, we know that ^(.^i?) will be taken as zero. 
Define 


hi(E’) = 

/ (T s , (E,EA) dE 

’Ei 

(Ej< E’ <E i+ 1) 

(67) 

hrAE’) = 

f E i+ 1 

/ *„(£,£') dE 

)u\E' 

( fi'/H-l < E’ < Ej+i) 

(68) 

h(AE') = 

f E> 

/ <7^ (£,£’') dE 

JEi 

(Ej < E’ < E l+ 1) 

(69) 

h <U)^ = . 

f Ej-\-j + 1 

/ <TsJE,E') dE 

J E, +J 

(A ( ; +J <^.<A (+;+1 ) 

(70) 

W*') = 

rE' 

/ a ri AE,E')dE 
JEj 


(71) 

*4(d)(^ = 

f E ' + ' (r ry i}(E,F?) dE 


(72) 
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then write the coefficients associated with the system of differential equations as 


a i,i — ^4 + ^6 + + ^3(2) ~ ^ ' 

a »\»+l = ^5 + *7(1) + *4(1) + *4(2) 
a M+ 2 = *7(2) + *4(1) + *4(2) [ 


( 73 ) 


where evaluation at the appropriate mean energies is implied. 

In this way, we generate a system of equations having the triangular form given by 
equation (34). We again use the source terms g(Ej,Xj ,y k .) obtained from the HZETRN simulation 
of the solar particle event of February 23, 1956, associated with an aluminum- water shield- 
target. configuration. Note that now we must solve the multigroup equation (34) associated with 
equation (40) for the multiple atom target material of water. We consider the cases of discrete 
shield thickness a?2,*3>— and apply the multigroup method to the solution of equation (16) 
applied to all target material y > 0. For each revalue considered, the initial conditions are 
obtained from the previous solutions generated where y — 0. This represents the application of 
the multigroup method to two different regions: region 1 of all shield material and region 2 of all 
target material. We then continue to apply the multigroup method to region 2 for each discrete 
value of shield thickness, where the initial conditions on the start of the second region represent 
exit conditions from the shield region 1. This provides for continuity of the solutions for the 
fluence between the two regions. 

Reaction Effects on Evaporated Neutron Fields 

In the present calculations, the two-stream bidirectional version of the multigroup method is 
always used because of its improved physical description and improved accuracy, especially near 
the boundaries of the incident radiation. Here, the assumption is made that half the evaporation 
source neutrons move in the forward direction and the other half move in the backward direction. 
The multigroup equations are, therefore, solved twice, once for the forward half of the source 
term and again for the backward half of the source. We evaluate the radiation fields for the 
solar particle event of February 23, 1956, in an aluminum slab 100 g/cm 2 deep with the results 
shown in figure 4 using the computational code of Heinbockel, Clowdsley, and Wilson (ref. 19) in 
which the evaporation neutrons are transported under elastic scattering only. Also shown in the 
figure are results obtained with the nonelastic process as described by the present calculation. 
A general decrease occurs in the 5 to 25 MeV neutron flux with a corresponding increase below 
2 MeV. As one would expect, the more reactive energetic neutrons are removed from the field 
by the reactions with the appearance of lower energy neutrons as reaction products. Also shown 
are results from the MCNPX Monte Carlo code. It is clear that the discrepancies reported by 
Clowdsley (ref. 22) are not from neglect of nonelastic processes. Similar results are shown in 
figure 5 for a water target along with Monte Carlo calculations using the LAHET code (ref. 6). 
The discrepancies observed in our earlier calculations are clearly due to factors other than effects 
of reactive processes associated with the transport of evaporation neutrons. 

High-Energy Backward Produced Neutrons 

Although the two-body interactions of nucleons are limited to the forward scattering, the 
multiple scattering of nucleons in nuclei can produce a nucleon in the backward direction after 
several scattering events. In addition the Fermi motion within the nucleus will enhance this 


16 



effect. The angular distribution of nucleons from nuclear reactions was estimated by Ranft 
(ref. 23) to be given by the approximate function 


g(A h Eft) = 


-0 Z 


N exp 


N exp / — — 


2 

"ET 


(o< e < tt/2) 

(Otherwise) 


(74) 


where A = (120 + 0 M6Aj^)/E, with E the secondary particle energy in MeV, Aj the atomic 
weight of the struck nucleus, and N a normalization constant. The fraction of neutrons produced 
in the forward direction is 

F[ or = 27t f g(ATiE,0) d{cos0) 

J for 

The corresponding backward-produced neutron fraction is 


(75) 


E bac = 1 - F[or 

The approximate isotropic component of the interaction can be taken as 

F\so — 2 F\ Y(ir 


76) 


77) 


In earlier development of the multigroup method, we assumed that the direct reaction products 
would mainly be of high energy and in the forward direction and, therefore, adequately sol ved 
by the HZETRN code. The assumed isotropic evaporation source was then treated with 
the multigroup method by assuming half of the evaporation source was propagating in the 
forward direction and the second half in’ the backward direction. In similar fashion, we replace 
evaporation and direct reaction spectra in the code as follows: 


f*(E,Ef) - F\*Jif'(E t E') + f d (E,F?)] 
f d (E,£f) — i'HZETRNLT (E,E') + f d {E,E’)] 


(78) 


where TjjzeTRN = 1 — F^ 0 . These replacements (eqs. (78)) were made in the new version 
of HZETRN/ multigroup code which is now only a minor modification. The terms on the 
right-hand side of equations (78) are shown in figure 6 and should be compared with figure 1. 
The importance of the reactive channels are accentuated because of the higher energies of the 
backward propagating neutrons. 

Results for Ranft Modified Source 


We have reevaluated the neutron fields in aluminum and water for the solar part icle event 
of February 23, 1956, with the angular dependence of Ranft and the separations into HZETRN 
and isotropic components. The results are shown in figures 7 and 8 along with MCNPX (ref. 22) 
and LAHET (ref. 6) derived Monte Carlo results. The addition of the high-energy backward 
component is essential in reaching agreement with the MCNPX and LAHET codes. It is clear 
that the discrepancies observed by Shinn e t. al. (ref. 16) in the 50 to 200 MeV region are due 
to the energetic neutrons produced in the backward direction as reasonably described by the 
Ranft formula. The Ranft formula appears to overestimate the backward component for oxygen 
because' agreement is improved for the omnidirectional flux at larger depths although forward 
and backward components may be somewhat incorrect . Quite satisfactory agreement is obtained 
at the largest depths. Still many of the cross sections in the HZETRN code are crude and a 
continued effort to improve them is expected to further enhance the calculated results. 
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Improved Cross Sections 


A program of improved cross sections has been in progress for several years. Greatest 
attention has been given to fragment events for heavy ions which have been significantly improved 
(refs. 24 and 25). Recent years of research have resulted in improved absorption cross sections 
(ref. 26) and improved production cross sections for the present study using the LAHET code 
with results in table 2. The effects of these new cross sections are shown in figure 9. 


Table 2. Number of Nucleons Produced in Nuclear Collisions With Aluminum Atoms 


Energy, MeV 

Cascade Nucleons 

Evaporation Nucleons 

n — - p 

n — n 

p — n 

P~*P 


n — ?? 

p — ^ n 

P^P 

25 

0.13 

0.26 

0.09 

0.24 

0.38 

ME EH 

0.43 

— 

200 

0.73 

1.26 

0.80 

1.18 

0.75 

wEm 

0.90 

IfM 

400 

0.91 

1.58 

0.99 

1.49 

0.87 

mm 

1.04 

mm 

1000 

1.48 

2.12 

1.59 

2.02 

1.30 

1.69 

1 .55 

1.44 

. 

1.97 

2.58 

2.11 

2.45 

1.44 

1.83 

1.70 

1.58 

1 

2.32 

2.92 

2.49 

2,76 

1.48 

1.86 

1.72 

1.62 


Concluding remarks 

The methods described herein greatly improve the HZETRN computer code's neutron 
transport predictions. To summarize, a bidirectional multigroup solution of the straight-ahead 
Boltzmann equation for elastic and nonelastic transport of low-energy evaporation neutrons has 
been implemented. The resulting computer code was added to the existing HZETRN computer 
code which was developed at the Langley Research Center. With the new modified code, various 
simulations were conducted to test its accuracy. The Monte Carlo codes LAHET and MCNPX 
were used as benchmarks of accuracy. The modified code with and without the inclusion of 
nonelastic scatt ering processes for the evaporation neutrons is compared with these benchmarks. 
The neutron fluences are calculated at depths in aluminum of 1, 10, and 30 g/cm 2 . These depths 
were selected because they represent typical values of shielding associated with the constantly 
changing space environment encountered by astronauts. The shield material of aluminum is 
typical because of weight considerations in space. The neutron fluences are calculated at 
depths of 1, 10, and 30 g/cm^ in water. Water is used to model human tissue. Including 
nonelastic scattering processes in the calculation of the transport of low-energy evaporation 
neutrons slightly improves the prediction of neutron fluence at low energies, but the prediction 
of neutron fluence at slightly higher energies, around 10 MeV, is decreased by this change. 

The HZETRN code underestimates the fluence of neutrons in the range of 5 to 200 MeV. In 
an effort to fix this problem, a formula by Ranft was used to estimate the number of isotropic 
neutrons at each energy. Using the bidirectional multigroup method to propagate all the isotropic 
neutrons greatly improved the neutron fluence predictions. 

The addition of improved production cross sections to the code showed only modest improve- 
ment to the predicted neutron fluence. In the future, more accurate absorption cross sections 
will be added to the code, but this is expected to also have only a modest effect. 
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Figure 1. Evaporation and direct cascading ne 
neutrons in aluminum. 
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Figure 3. Multigroup energy partition. 


23 




10' 2 I0' 1 10° 10 1 10 2 lO 4 10 4 

Energy, MeV 


(a) Depth of 1 g/crn^. 



Energy, MeV 

(b) Depth of 10 g/cm 2 . 

Figure 4 . Energy spectra of neutron fluence in aluminum calculated by HZETRN program with 
bidirectional multigroup method used to transport evaporation neutrons. 
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Figure 6. Forward moving and isotropic neutron spectral effects for collision of 500 MeV neutrons 
in aluminum. 
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Figure 7. Concluded. 
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